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ABSTRACT 


A numerical study was performed to investigate the dynamic 
crack propagation in fibrous composite plates utilizing the 
finite element method. A rectangular plate of uniform 
thickness, which had a propagating central crack, was used for 
the study. The plate was a unidirectional composite panel and 
the load was applied in the longitudinal direction of the 
composite plate. Fracture energies were calculated for given 
speeds of cracks. The objective of the study was to examine 
the effect of different composite material properties, crack 
speeds, and densities on the fracture energy. The 
discontinuous node-release technique was used to model the 
crack propagation. 

The numerical study showed that the fracture energy was 
higher for a lower elastic modulus ratio, if all other 
conditions were held the same. Furthermore, a lower crack 
propagation velocity or a lower material density resulted in 
a higher fracture energy, respectively, provided the rest of 


the parameters held constant. 
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ibe INTRODUCTION 


Dynamic fracture deals with problems having a crack 
within a body where the inertial forces play an important 
role. The early investigation of the dynamic fracture analysis 
goes back to the case of unstable crack propagations in some 
ships during World War II. The dynamic fracture problems 
involve crack initiation and crack propagation. There are two 
principal types of problems. One is applying a dynamic load to 
a body containing stationary cracks, and the other is 
concerned about a body with a moving crack. The second type of 
problems can be classified into two different modes; 
propagation and generation modes. 

Stationary crack problems involve cracks within structures 
that are subjected to dynamic loading, such as an impact load. 
For example, Chen performed a numerical analysis to find the 
dynamic stress intensity factor for a rectangular plate having 
a central crack. He used the finite difference technique and 
discussed the dynamic behavior of the stress intensity factor. 
He argued that the oscillations on the stress intensity versus 
time curve were mainly due to scattering phenomena from the 
crack tip and the boundary surfaces. He concluded that the 
overshooting of the dynamic stress intensity factor is 


attributed mainly to the geometry and loading. [Ref. 1] 


The generation mode involves moving cracks expanding at 
prescribed rates. The main objective of the study is to 
compute the stress intensity factor or the fracture energy and 
to determine the fracture toughness. The propagation mode 
deals with the cases where fracture toughness is known. The 
objective is to predict the crack propagation and its arrest 
for a given load. One of useful models to investigate the 
crack propagation and arrest was the double cantilever beam. 
Kanninen and his coworkers [Ref. 2] proposed the model and 
Studied it extensively. Some of studies in the dynamic 
fracture were listed in the references [Ref. 3-7]. 

Although there are many literatures available for the 
dynamic fracture of isotropic materials, the studies of 
composite materials are not so common. Some of the fracture 
studies of composites were shown in references [Ref.8-10]. 

This study considered a rectangular plate made of 
unidirectional composite and containing a central, moving 
crack subjected to a uniform tensile load. The finite element 
method utilizing bilinear rectangular elements was used for 
this study. A parametric study to compare fracture energies in 
composite plates for various conditions was conducted. 
Different elastic modulus ratios, crack velocities, and 


material densities were considered. 


Il. MATHEMATICAL DERIVATION 


A. GOVERNING EQUATIONS 

The governing equations for the problem of elasticity are 
equations of equilibrium, constitutive equations, and 
kinematic equations. For the two-dimensional problem, they are 


given below. Equations of equilibrium are 


Beret, fu 
Ox dy ~ gt? 
iy RO Ga 
Ox oy at2 








(1) 








and T are the two normal stresses in the x and 


where 0,,0y, xy 


y direction and the shear stress, respectively. p is the 
material density, and u and v are the displacements in the x 
and y direction, respectively. In addition, x and y are 
rectangular coordinate axes and t denotes time. 

Constitutive equations are also known as stress-strain 


relations. The relations are given in a matrix form like 


{o}= [D] {e} (2) 


where {co} and {e} are, respectively, stress and strain 


vectors, that is 


{o}=(0, o, t,J7 
lel=le, €, Y,,j7 


(3) 


and [{D] is the material property matrix of size 3x3. For an 
isotropic material, the [D] matrices for plane stress and 


plane strain conditions, respectively, are 








1 vu 0 
Do ee (4) 
(LSD ls, (lee 
2 
and 
it e 0 
aU 
= E(1-v) 1») ah 0 
[D] (1 eee (5) 
2 U 
| : ° 241 =U» 


For a composite material, [D] matrix becomes 


Di, Diz Dy; 
[D] =|D2, Dz2 Do (6) 
D3, D3, Ds; 


The details of Eq.(6) are shown in references. [Ref. 11] 
Kinematic equations are also known as strain-displacement 
relations. The infinitesimal strains are expressed in terms 


of displacements as shown below: 


Ou 


Ox 
ex 
“}- = (7) 
Na) | du, dv 
Oy Ox 


B. BOUNDARY CONDITIONS 

On the essential boundaries displacements u and/or v are 
prescribed, and on the nonessential or natural boundaries 
tractions are prescribed. Usually on a part of the boundary 
tractions are prescribed, while on the remainder of the 


boundary displacements are known. The tractions are given by 


1 ecole Wem = 8 
Da yt Oy 


(8) 


where 9%, and ey, are the tractions in x and y directions, 
respectively, and n, and ny are, respectively, the components 
of outward unit normal vector to the boundary. 

Applying the method of weighted residual technique to 


the equations of equilibrium yields 


Ou Oo, 


8x, iy 
8t2 ox oy 
7 Gy ot Oa . 
a eT Des ae aa 


T,=Jg(-p 








) w,dQ=0 
(9) 





where ow, and w, are test functions and 2 is the domain of a 


given problem. Applying the weak formulation results in 


Ou dW, JW, 
mee Fl, +1 
= 8 es Y dQai yan = 9 (10) 
p ——_ W,, +t, ~~ +9,-~— 
$e" YT OX ae 


and fT is the boundary of domain ©. Rewriting Eq.(10) gives 
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are the mass matrix, stiffness matrix, and load vector, 


respectively. 


A finite element discretization is applied to the domain 


and the following matrices are computed at the element domain. 


C. MASS MATRIX 


The mass matrix is given by 





Ome. 
0, 0])9e2 

1 

Jar P| g be 32 ,,( AR - (15) 
Ce 


In this study both bilinear rectangular and linear triangular 
elements will be used to produce the finite element mesh. 
First the finite element derivation will be shown in detail 
using the bilinear rectangular element. The rectangular 


element is shown in Figure 1. 





Figure 1 Bilinear Rectangular Element 


The shape functions are 








See Z 7 
H, => 22 x) (c-y) 
= 1 — 
BG TEE Ot) es) han 
1 
H,= nyse (b+ x) (Eig) 


ee 2 
Hy = Fae (D-x) (c+y) 


where b and c are greater than zero. The displacements in the 
x and y directions, respectively are interpolated as shown 
below 

u=H, U, +H,U,+H,u,+f, u, 


(17) 
V=H) Vato Vo ave, VV, 


where 


2 (18) 


is the displacement vector at the nodes of the element. 
Because the shape functions are independent of time the 


acceleration terms become as shown below: 


07 u 
eo) af uy +H, U. +H, U, el 
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op (19) 
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Subst i Elem ekope@mc back anto Eq.(15) and applying 


Galerkin’s method results in 


H, 0 
0H, uy 
H, 0 a 
U2 
ame 8 8 OF Ce, ony 
1H, OF0 H, O H, O H, 0 H,! 4s 
V3 
0 4H, i 
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Finally Eq. (20) becomes 


Hi (0 Hj, 0° 4,4, (ORE 0 
0 Hi 0 TH 0° nie HOF 


Uy 
H,H, O Hy; OO 4H,H, O 4H,H, 0 Vy 
Up 
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fase x i Naat (22) 
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H,H, O HH, O 4H,H, O Hy 0 
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Carrying out the integrations gives symmetric consistent 


mass matrix as shown below: 


~~ © 
mm OC DW 
Pf O NY © 


(22) 


Cc 
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[M.] = 2 be 


Pm ON OF O 
PrP ONM OF OO ND 
PON OF ON O 


where t is the thickness of the element. Equation (22) is 
called the consistent mass matrix. 

Another choice for the mass matrix is the lumped mass 
matrix. The lumped mass matrix is a diagonal matrix so that it 


reduces computation time and space needed for storage. The 
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lumped mass matrix is 











iO..0) 0 0: 0: 0 C0 
moO 0 0 0 0 
10000 0 
1 070)-0 0 
M,)=pbc 
[M,] =p 1000 
iOr0 
jb 18, 
1 
D. STIFFNESS MATRIX 
The stiffness matrix is given by 
dW, OW, : 
Ce  . 0 _ 
le * oda . 
dW, dW,|I, 
Oy Ox 
Substituting constitutive equations, Eq. 


yields 
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(2), 


0 | 


xy 


(23) 


(24) 


ma cO sg .24) 


(25) 


Using kinematic equations, Eq. (25) becomes 
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Xx Y OV 
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The kinematic equation can be written in terms of the nodal 


displacements as below: 





du 
Ox oO 0 ne) 
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where the last matrix is called [B] matrix and the last vector 
called the displacement vector {u}. Substituting back Eq. (27) 


into Eq. (25) yields 











OW, OW, 
Ox oy 
- : (28) 
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Oy Ox 


Using Galerkin’s method the first matrix in Eq. (28) becomes 
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ag at 

Ox Oy 
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Ox Oy 
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which is equal to the transpose of B matrix. 
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Finally, the stiffness matrix is 


[Kl =[o- [B]7(D] [Bl dQ. (30) 


Computing the [K] matrix using bilinear rectangular element 


gives 


[5.f?,(B] *(D] [B] tdxdy (31) 
where t is the thickness, and 
=e, 0 (c-y) 0 (ct+y) 0 -(c+y) 0 
[Bl =| 0  -(b-x) 0 -(b+x) 0 (b+tx) 0 (b-x)| (9) 
=({b-x) =(e-y) =(ee)) (e=y) pibtx) (cry) (bac ecay, 
For a general material property matrix [D] of Eq.(6), the 
element stiffness matrix becomes 
Ky, Ky Ky 3 Ki, Ky 5 Ki Ky Ki 
Koy Ky K33 Ko. Kos Ky kK, Ko 
R31 K3> K33 K3, K35 K36 K37 K3 3 
[K] Kyy Ky Ky, Ky, Ky45 Ky¢ Ky Ky , 
a 33 
Ks1 Ks Ks, Ks4 Ks- Ks Ks7 Ks 
K; al Ke i Ke 3 Kea Ke 5 Ke 6 K67 Ke 8 
Koy Koo kK, Ko, Ks Ko Koy og 
| Kg ah Kg a Kg ie} Kg 4 Kg 5 Kg Kg7 Kg 8 


The detail of K,;; is shown in the Appendix. 
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If a material is isotropic and the plane strain condition 


applies, the stiffness matrix is 


[K] =(K,]+([K,] . (34) 


[K,] is given by 


= bic 4ab*c? = bic -4dab*c? bc -4ab*c? "bic 4ab*c? 

=< be? aab2c? —48bc3 -4ab2c? —8bc3? -4ab 2c? $8 be3 
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16503 gab2c2 218 be3 

3 5 
+2 bc -Aab*c? 
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where 


. Eyl) oa 
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(36) 
v 
(=) 


— 
= 


and [K,] is given by 


als, 


+2 bic Ab*c? Sb3c “4p -¢" bic -4b*c? bic 4b*c? 


16 be3 ab2q2? 18 pbc} -ab2c? Bbe3 -ab2c? = 8be3 
3 3 3 3 
= bic -4b*c? =" b3c -Ab*c? bic 4b*c? 
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3 3 5} 
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<2 be? 4b2c2 ==" be? 
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= be? 
where 
p= ee (38) 
Pa oe) — (lly) 
E. LOAD VECTOR 
The load vector is given by 
“iy 
I . (39) 
ig 


Linear shape functions can be used to find nodal load values 


due to the traction. This results in lumping tractions into 
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two nodal points associated with the boundary of the element. 
ihtmehere 19)no —traceton, tue lead vector is zero at the 


associated nodal points. Linear shape functions are 





Seques 
arn 
oh 
SS) (40) 
H, = a 


aa 
As=S4447S; 


where s is the coordinate axis along the boundary, and h,; is 
the length of the boundary of the element. Combining all of 


the above procedure gives an equation of matrix, 


[M] {t+ [K] {ub={R} (41) 


which describes the equation of motion with no damping. 


Damping was neglected in Eq. (41). 


F. LINEAR TRIANGULAR ELEMENT 

Since this study deals also with the linear triangular 
element, the derivation of the matrices will be shown here, 
but not in detail since the detail was shown before in the 
bilinear rectangular element case. The linear triangular 
element, as Shown in Figure 2, has three nodal points which 
have coordinate values (X,Y), (X2,Y2), and (X3,Y3) 
respectively. Linear triangular element has two degrees of 


freedom per node. 


iL 





Figure 2 Linear triangular element 


Note that three nodal points are numbered in the counter- 
clockwise direction from an arbitrarily selected corner. 


The shape functions are 


i 
A, = 3A [A, +B, xX+ yl 
a [A, +B,x+C,y] (42) 


= [A, +B, x+C,y] 


where A is the area of the triangle, in other words 


1 x, Vi 
A=s det j1 x, y2 (43) 
ey Oe 


and 


dite: 


A, a 2am 3) 2 B, myo ays ei =X, — Xo 
Az =X3¥1~*X1V3 B.=¥3-Vi CSS 8S (44) 
A, =X, V2~-X2Y, By=¥1-y2 ee ia ome 


It can be noted that 
A=(A,+A,+A,)/2 . (45) 
The displacements in the x and y directions are given by 
USE, HU, FH, (46) 
Wet Vaid V5 id, Va 


Carrying out similar calculations and procedure gives the 


symmetric consistent, and lumped mass matrices, respectively 


(47) 
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The stiffness matrix is obtained from a combination of two 


components 


ko 
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2 
oe EV Ci Gos, 6, C,, Go: (50) 
* 8A2(1+v) Be B,C, Bees 
Cl GBs 
3 
where V is the volume of the element, that is, 


and t is the thickness of the element. 


G. TRANSIENT ANALYSIS 


In this study the central difference method is used as a 


time integration scheme. The equation of motion at time t is 
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[M] {ti} + (C]{u)*+(k] {udt=tr}e (52) 
The acceleration at time t is given by 
ti"-=—— SNe Alujt+{u}eree] (53) 


The velocity at time t is 


ere ally tb) label 0 handel fa (54) 


The error associated with the acceleration and velocity 
Goleulationasareseneerder (At-). Substituting the relations 


for acceleration and velocity into the equation of motion, 


Eth GC — 
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The method requires the solution at time -at. Using the 
equations of acceleration and velocity, letting time to be 


Paula ktouzenouaid eliminating {u}** yield 
2 
(ule =(u}—9 tled?+ 22) (ajo (56) 


Note that {ti}° can be found using the equation of motion, Eq. 
(52), at time zero since {u}° and {u}° are known. 
The central difference scheme is conditionally stable and 


requires that the time step at should be smaller than a 


Pgh 


critical time step at... The critical timegsreo seme ivory, 


ati.= = (57) 





where T,., is the smallest period. 


n 


Therefore the method can be summarized as; 
1.Compute mass, damping, and stiffness matrices. 


2.Compute {a}9 from 
[mM] {tx}? ={R}°- [ C] {uz} + [K] {ul}? (58) 


where {u}-,and.4u} are Umertailconcus tomer 


3.Compute 
{y}-at={1y}0-, t{u}e+ 1A£)” (ajo | (59) 


4.Compute effective mass matrix from, 


— ——— pe 
(A = IM + cl (60) 


5.Compute the effective load vector at time t, 








ta{phe_ = eee t_ a oo ae t-at 
{Rif ={R}- ( [X] ey 2 le ‘Feng a Lela See) 


6.Solve for displacements at time t+at from, 


[M) {ubttee={R}e (62) 


7.Compute accelerations and velocities at time t, 


Z2 


(pea tm Pea {ule (ules) 
a (63) 
= ae ene eee 
{ ae ] 
Note that steps 5-7 have to be performed repeatedly for each 


time increment. 


H. DYNAMIC FRACTURE ANALYSIS 

There are many parameters to describe the crack tip 
behavior for determination of the crack propagation and 
arrest. One of the parameters is the stress intensity factor, 


K, at a crack tip which is defined as 


Rowen 2G (ten 0r Gc) (64) 


where 


CP =OmCanU, C) (65) 


is the stress component, and is evaluated near the crack tip. 
In this case, crack propagation is governed by stress 
Hisenstiay racu@rem caiduehe material property @K, , calledwethe 
fracture toughness. A crack becomes unstable if the stress 
intensity factor reaches the fracture toughness of the 


material. 


ZS 


The dynamic analysis requires K be considered in two 


parts. For crack antetoniom 


K(a,o, €) =Kp(o) (66) 


where a is the crack length, o is the applied stress, t is 
time, and o6 represents the loading rate. For a propagating 


crack, 


K(a,o, t) =K,(4a) (67) 


where a indicates the crack speed. [Ref. 7] 


An alternative criterion for crack motion is, 


GilenG bs) h(a) (68) 


where G is the dynamic energy release rate and R is the energy 
dissipation rate required for crack growth which is the 
Critical value of G. 

The driving force or the dynamic energy release rate for 
crack propagation has contributions from strain energy, 
kinetic energy, and work done on the body. The driving force 
is the net change in these components per unit area of crack 


extension. The equation for G is, 
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aW duU_dadT 


G= jl A ald alia 


$ da da da 


“44 (dw_du_ ar, ee 
“ba dt dt dt 


where U is the strain energy, T is the kinetic energy, W is 
the work done by the external loading, a is the crack length, 
and b is the plate thickness at the crack tip. 

The dynamic energy release rate can be directly connected 
to the dynamic stress intensity factor. For plane strain 


conditions, this relation is 


g- 12% a(4) K (70) 
1h 


where E and v, as usual, are the elastic modulus and Poisson’s 
Paedo,e reSpeceiveiy mand) A 1S a geometry independent function 


of the crack speed given by 


a 2 Se he 3 
iofeed? (1-3) 
A= : 
(1-v) [4(1- 42) 7 (1-20) F- (2-22) 2] 
Or Cy C3 
where C, and C, , respectively, are longitudinal and shear 


(transverse) wave speeds in the material. Wave speeds are 


given by, 


1/2 
C,= aa (72) 


p 





Zi 


and 


@=(G/p (73) 


where K, and G, are bulk and shear modulus, respectively: 
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IIil. RESULTS AND DISCUSSION 


A. VERIFICATION PROBLEMS 

In order to verify the developed computer code, two 
example problems with known solutions were analyzed for 
isotropic plates. One was a stationary crack problem subjected 
to a rapidly varying dynamic load, and other was a propagating 
crack problem with a constant speed and subjected to a uniform 
load. 

The first verification problem was a stationary crack 
problem known as Chen’s problem [Ref.1]. A rectangular bar 
with a centrally located crack, shown in Figure 3, was 
considered. A uniform tension of magnitude 0.004 Mbar with 
Heaviside-function time dependence was applied. The material 
properties were given as below; shear modulus of 0.76923 Mbar, 
bulk modulus of 1.66667 Mbar, Poisson’s ratio of 0.3, and 
material density of 5 gr/cm*®. Chen solved the problem using 
the finite difference technique. 

Both Chen’s solution and the present solution are shown in 
Figure 4. A good agreement between the two solutions was 
obtained. 

As the second example a rectangular plate with a central, 


moving crack was considered [Ref. 3]. The crack was assumed to 
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Figure 3 The Geometry of the Stationary Crack Problem 
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Figure 4 The Results of the Stationary Crack Problem 
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propagate at a constant speed with uniform stress field of 
68.9 MPa. The crack speed was assumed to be 33% of the 
dilatational wave speed in a steel which is 5904 m/s. The 
material properties were elastic modulus of 206.7 GPa, 
Poisson’s ratio of 0.30, and material density of 8000 kg/m?. 
The geometry of the problem is shown in Figure 5.([Ref. 3] 

A good agreement was obtained between the analytical and 
the present solutions. Both the analytical and present 


numerical analysis results are shown in Figure 6. 


B. RESULTS OF COMPOSITE PLATES 

A composite plate with a central, moving crack, with a 
Similar configuration as the previous example, was 
investigated for different cases. A unidirectional composite 
plate was considered. 

For all cases calculations were made for work done on the 
system, internal energy, kinetic energy and fracture energy. 
The objective of this study was to compare the fracture 
energies for various conditions. The first case was the 
variation of the ratio of elastic moduli of the longitudinal 
and the transverse directions while keeping the others fixed. 
The next case changed crack velocities and, the last case 
varied material densities while keeping the other parameters 


fixed. The geometric and material properties of the problem 


2G 


wre ‘7 ra 


i | ||| 
eee eee 
Pett eet tty 
PET Pf 12m 
Pet ETT yt 


= 4 aes 


Crack tip 





Figure 5 Finite Element Mesh for the Moving Crack Problem 
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Figure 6 The Results of the Moving Crack Problem 
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are length of 3m, width of Im, elastic modulus ratio in the 
transverse direction of 10 GPa, Poisson’s ratio in the 
transverse direction associated with loading in longitudinal 
direction of 0.25, and uniform loading in the longitudinal 
Ginectien of 70 MPa® 

Figure 7 shows work and energy variations during a crack 
propagation. For all cases throughout this study, work was 
higher than internal and kinetic energies, and internal energy 
was higher than kinetic energy. 

In the first case, the elastic modulus ratios were 
changed. Three different ratios between the elastic moduli in 
the longitudinal and transverse directions were used. The 
crack velocity and the material density were kept constant. 

Figure 8 and Figure 9 show work and internal energy 
variation, respectively. The results reflect higher work and 
higher internal energy for a lower elastic modulus ratio. A 
softer material which has a lower elastic modulus ratio is 
subjected to a larger displacement under the same load while 
all other parameters held the same. This produces higher work 
and higher internal energy since they are proportional to 
displacements. As will be seen in all cases, the trends of 
work and internal energy are always similar to each other. 

Figure 10 shows the kinetic energy variation. Kinetic 


energy is higher for a lower elastic modulus ratio, because 


ee 


for this material the material velocity is higher due to the 
larger displacement. Figure 11 shows the fracture energy 
variation. Fracture energy is also higher for a lower elastic 
modulus ratio. 

The increments of work and energies between the initial 
time and a later time are plotted in Figure 12. The larger 
fracture energy for a lower elastic modulus ratio was caused 
by a much larger increase of work relative to the increases of 
internal and kinetic energies for this material. 

In the second case the crack velocity was varied. Three 
different crack velocities were used. The elastic modulus 
ratio and the material density were kept constant. Figure 13 
and Figure 14 show work and internal energy variations, 
respectively. The time needed for a crack propagation to the 
same crack size increment is longer for a lower crack 
velocity. With a longer time, the crack propagation influences 
the rest of the body more widely. The wider influence results 
in a larger displacement. As a result, it produces higher work 
and higher internal energy. A typical displacement variation 
at a point, shown in Figure 15, clearly is larger for a lower 
crack velocity. 

Figure 16 shows the kinetic energy variation. Kinetic 
energy is higher for a higher crack velocity. The reason is, 


as shown in Figure 17, a higher crack velocity produces higher 
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speeds in the material. Figure 18 shows the fracture energy 
variation. Fracture energy is higher for a lower crack 
velocity because a lower crack velocity caused a larger 
increase of work and a less increase of kinetic energy than a 
higher crack velocity. 

The material density was varied as the last study. Three 
different material densities were used. The elastic modulus 
ratio and the crack velocity were kept constant. Figure 13 and 
Figure 14 show the work and internal energy variations, 
respectively. For the same magnitude of an external force, a 
smaller density causes a higher velocity and larger 
displacement of the plate. Because the displacement is larger 
for a lower material density, work and internal energy are 
Baaher:. Figure 22 confirms it. Figure 23 also shows a higher 
speed for a lower material density. 

Figure 24 shows the kinetic energy variation. Kinetic 
energy is higher for a higher material density. The velocity 
is higher for a lower density, but the mass is larger for a 
higher density. Because the kinetic energy is proportional to 
the mass as well as the square of velocity, which case gives 
a higher kinetic energy depends on which term is dominating. 
In this situation large mass yields a higher kinetic energy 
although the material velocity is lower. Figure 25 shows the 


fracture energy variation. Fracture energy is higher for a 


so 


lower material density because of the same reason as that for 
the case of different crack velocities. The increments of work 


and energies are also plotted in Figure 26. 
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Figure 7 Work and Energy Variation During a Crack Propagation 
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Figure 8 Work Variation for Different Elastic Modulus Ratios 
and Fixed Crack Velocity and Material Density 
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Figure 9 Internal Energy Variation for Different Elastic 
Modulus Ratios and Fixed Crack Velocity and Material Density 
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Figure 10 Kinetic Energy Variation for Different Elastic 
Modulus Ratios and Fixed Crack Velocity and Material Density 
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Figure 11 Fracture Energy Variation for Different Elastic 
Modulus Ratios and Fixed Crack Velocity and Material Density 
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Figure 12 Work and Energy Differences Between the Initial Time 
and a Late Time for Different Elastic Modulus Ratios and Fixed 
Crack Velocity and Material Density 
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Figure 13 Work Variation for Different Crack Velocities and 
Fixed Elastic Modulus Ratio and Material Density 
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Figure 14 Internal Energy Variation for Different Crack 
Velocities and Fixed Elastic Modulus Ratio and Material 
Density 
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Figure 15 Displacement Values for Different Crack Velocities 
and Fixed Elastic Modulus Ratio and Material Density 
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Figure 16 Kinetic Energy Variation for Different Crack 
Velocities and Fixed Elastic Modulus Ratio and Material 
Density 
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Figure 17 Speed Values for Different Crack Velocities and 
Fixed Elastic Modulus Ratio and Material Density 
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Figure 18 Fracture Energy Variation for Different Crack 
Velocities and Fixed Elastic Modulus Ratio and Material 
Density 
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Figure 19 Work and Energy Differences Between the Initial Time 
and a Late Time for Different Crack Velocities and Fixed 
Elastic Modulus Ratio and Material Density 
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Figure 20 Work Variation for Different Material Densities and 
Fixed Elastic Modulus Ratio and Crack Velocity 
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Figure 21 Internal Energy Variation for Different Material 
Densities and Fixed Elastic Modulus Ratio and Crack Velocity 
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Figure 22 Displacement Values for Different Material Densities 
and Fixed Elastic Modulus Ratio and Crack Velocity 
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Figure 23 Speed Values for Different Material Densities and 
Fixed Elastic Modulus Ratio and Crack Velocity 
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Figure 24 Kinetic Energy Variation for Different Material 
Densities and Fixed Elastic Modulus Ratio and Crack Velocity 
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Figure 25 Fracture Energy Variation for Different Material 
Densities and Fixed Elastic Modulus Ratio and Crack Velocity 
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Figure 26 Work and Energy Differences Between the Initial Time 
and a Late Time for Different Material Densities and Fixed 
Elastic Modulus Ratio and Crack Velocity 
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IV. CONCLUSIONS AND RECOMMENDATIONS 


The developed computer code has been proven to work well 
for dynamic crack problems. Very good agreements were obtained 
for the known solutions. The parametric study of dynamic crack 
propagation in composite plates provided the following 
conclusions. 

A softer composite plate which has a lower elastic modulus 
ratio showed a higher fracture energy if the crack velocity, 
the material density and other characteristics of the problem 
were held constant. A Benmos ite plate with a lower crack 
velocity and a composite plate with a lower density showed a 
higher fracture energy provided all other parameters were held 
constant, respectively. 

In this study a discontinuous crack release technique was 
used to model the crack propagation. This causes sudden shock 
waves which in turn result in oscillations in solutions. A 
technique that releases the propagating crack nodes 
continuously, for example as shown in Ref. 12, is recommended 


to achieve more stable solutions. 
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APPENDIX 


The general stiffness matrix using a bilinear rectangular 
element and the general case of material property [D] matrix, 
that is, 

Diy 
[D] =|P21 

D, 


1 


iS os 


N 


is given as below: 


Kyy Ky2 Ky; Kyi, Kis Ky¢ Ky7 Ki 
Koy Kp» Kp; Koa K5 R26 K27 Kz 
K31 K32 K3; K3, K35 K3¢ K37 K3¢ 
[K] " Kay Ky. Ky; Ky, Kgs R46 Kay Kas 
Ks4 Ks Ks3 Ks, Ks5 K5¢ Ks, Ks 
Kei Ke 2 Ke 2 Ke “hk Ke5 Ke 6 Key Ke 8 
Ky K K3 Rog Kas Ko Ky Kg 
Ke, Kg2 Kg; Kg, Kgs Kee Key Kea. 


where 


K,,==—D,,bo? +4 (Dp) b?c?+ = D,,b3c 


Ke ==2D, be? eB Balls) b?c?#+ =" D,,b3c 
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K,3#=-=2D,,bo? +4 (-D3,+D,3) b?c?+=D,,b3c 


K,g2- =D, bo? +4 (DD) b?c?+=D,,b7c 


K,s=-=D,,bo?+4 GD D..) b4 2—-—Dysb*c 
Kyg2-=D,,bo? +4 ( So E.), b*c?-=D,,b?c 
K.=2D..bc3+4 (D,.-D,.) b2c?-18p, bic 
17 a lal 31 13 oe 33 
K..=2D.,bo3+4 (-D,,+D,.) b2ce2-1© p,.b3 
18313 Cc 12733 Cc a wee Cc 


K,,==2 Dy, bo? +4 (D,,+D33) b?c?+="D,,b*c 


61 


Ky3=- =D, bo? +4 (-D,,+D,,) b?c?+=D,,b3c 


Kyg2- == Dy3be? +4 (De Dr) b?c?+=D,,b3c 


Kyg~ = D3,bC7+4 (—Dy~Dy3) b? 2- = Dy sb*C 


Kyg= 7 Dy,b07 +4 oO...) b?c?-=D,,b?c 


8 16 
Ky,=3D3,be°*+4 (oa D--)) b*c*-——D,,b*c 


8 iG 
Kyg=3 ya DC" +4 (ane) b*c*-——D,,b°c 


Ky,2- =D, be? +4 (DED) b2c#+=D,,b>c 


16 8 
Ky,==Dy,be*+4 (S10) bs) b*c*+—D;,b*c 


Ky" D,, bo? +4 (0b 0.) b?c?+ =" D,,b3c 
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16 


K,,=+°D,,bc?+4 (-D,,-D,,) b?c?++8 = D3,bc 


34 3 


Kyg==D,,b07 +4 (~D,,+D,3) b2c?- = D,,b°C 
Ky, = =D,; bce?+4 (D,,-D;,) b?c?- = Ds, bic 


Kyy=--=D,,bc7 +4 (DD, .) b?c?-2D,,b*c 


Kyg=-=D,,be7+4 coed.) b?c?-=D,,b*c 


Oe = = = Mba 4D s5= D33) b?c?+=D,,b3c 
K,,=-2°©D,,bc? +4 (-D,,+D,,) b?c?+= D,,b? 
a alee ad 320 220m eee 22 


Kar 310° +4 (-D,,-D,,) b*c* + D,,b*c 
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Kyg2 =~ D3,b0? +4 (2:50) b?c?+=" D,,b>c 


Kys==D3,b0? +4 (2D De) b?c?-="D,,b*c 


8 16 
Kye=Z D33be°+4 (Dae Des) b*c*-——D,,b*c 


8 


Ky,=- D3,be*+4 25+ D:-) b?c?-=D,,b°c 


Kyg2- = D3,b0? +4 (DD) b?c?-=D,,b>c 


K,,2--=D,,be? +4 (18 18) b?c?-=D,,bc 


K,g=-=D,,be3+4 (S10) 1B) b?c?-2D;,b°c 


Ko.=— Depo” AD) b?c?-=2D,,b3c 


woo 
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8 16 
Ky4=—D,,bc°+4 (GD o4 by,) b*c*-—=D,,b°c 


16 alte} 
K5=—-D,,bo* +4 (DED) b*c*+——D,,b°c 
6 ES 


ifs 16 
Ke oe el (10h seh, b*c*+—-D3,b°C 


16 | 8 
Kg=-—-D,,bc°+4 ELE Sek b*c*+—D3,b°c 


8 8 
K,.=-ZDybe* +4 (-D,,~D33) b*c*-—= D,,b°c 


Keg®~ = D3,bo? +4 (CD -De)) b?c?-=D,,b*c 


_ 8 16 
Ke3=3D3,be°+4 UDA, b*c*-——D,,b°c 
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8 16 
Keq=— D33bc° +4 (-D3,+D,3) b*c*-—— = DDE C 


— Ue 2,2, 16 
Kes=— a DC > Has aD EDC += Da3b*¢ 


Keg? == D,3be? +4 (Ds) Dee’ a Dep Ic 


a= = Dy, bc +a -Dz,+D,3) b?c#+=2 D,,b3c 


Keg=- = D,,bc2+4 (D,,—Dz3) b?c?+ =D, bc 

Kan ==D,, be3+4 (-D,,+D,,) b?c?- =2.D,,bc 
= 8p .bc3+4 (D,,-D,,) b2c?-+© p, pb? 

Spy Png Oe ee aces ori 


K,,2-=D,,bo7+4 (De coe b?c?-3D,,b3c 
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K,q2-=D,3be° +4 (D,,+D3,3) b?c?-=D,,b3c 


og ~ = —D,,bc*+4 (D,,-D,,) b*c? +g Dasb?e 


kag = +8 p,,bo3+4 (- -D,2+D33) b*c* +3Dsg oe 


Ky7= =D, bo? +4 (-D3,—D,3) b2c?+D,,b>c 


Kyg2= =~ D,,be* +4 ( -D, -D33) b*c*+—— - erate 


Kg, == D3,be3 +4 (-D,,+D;,) b*c*-—— = hae 


IG =8p, a oferta <0 Ree D,3) b?c?-=2 D,,b*c 


wl 


Kqy=-=D be*+4 (D,,+D33) b2c?-5 D,xb?C 


er 


8 8 
Kyq=-D33bc°+4 (D3,+D23) b*c*-— D,,b*c 


Kgs=- 2D, bo? +4 (D.. =D) b?c?+=D,,b%c 


Keg=- => Dy bc? +4 oe b?c?+=D,,b*¢ 


Kay “= , DC? Ha eo) b?c?+=D,,b°c 


Kyg2 =~ D3,bc? +4 (oo) b?c?+="D,,b3c 


In this study the [D] matrix for a unidirectional composite 


material for the plane stress conditions is; 


Di, Dy 9 
[D] =|D,, D,, 9 
0 O D,, 
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and the nonzero elements of the [D] matrix are; 


Ey 
ao 
1-U,,U5, 
EB, 
OO 
1-V,5V5; 
VizF,  _ U2 By 


BD =p ao a 
2 Zi = = 
iL U12V051 1 V1205, 


D,3=G,, 


where subscripts 1 and 2 on the right hand side of the 
equations are the direction of the fibers and the transverse 
direction, respectively. Substituting these values into the 
general stiffness matrix gives the stiffness matrix for a 


unidirectional composite. 
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